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We analyze the depletion interaction between two hard colloids in a hard-sphere solvent and 
pay special attention to the limit of large size ratio between colloids and solvent particles which 
is governed by the well-known Derjaguin approximation. For separations between the colloids of 
less than the diameter of the solvent particles (defining the depletion region), the solvent structure 
between the colloids can be analyzed in terms of an effective two-dimensional gas. Thereby we find 
that the Derjaguin limit is approached more slowly than previously thought. This analysis is in 
good agreement with simulation data which are available for a moderate size ratio of 10. Small 
discrepancies to results from density functional theory (DFT) at this size ratio become amplified 
for larger size ratios. Therefore we have improved upon previous DFT techniques by imposing test 
particle consistency which connects DFT to integral equations. However, the improved results show 
no convergence towards the Derjaguin limit and thus we conclude that this implementation of DFT 
together with previous ones which rely on test particle insertion become unreliable in predicting the 
force in the depletion region for size ratios larger than 10. 

PACS numbers: 05.20.Jj, 83.80.Hj 

I. INTRODUCTION AND PHYSICAL PROBLEM 

In many circumstances the effective interaction between larger colloidal particles suspended in a bath of smaller 
particles is determined by entropic effects. If both colloidal and solvent particle are taken to be hard spheres with 
radii i?2 and R\ = a/2, respectively, the colloidal interaction is purely entropic and arises mainly through the effect 
of a depletion zone between the colloids (for surface-to-surface minimal distances z < a) which is inaccessible to the 
solvent particles. 

For large size ratios a — R2/R1, the force F a in the depletion zone z < a between the colloids can be obtained 
by just using bulk and surface thermodynamics. This is the Derjaguin approximation (its derivation is presented in 
more detail below) which states: 

\ =P(*-g)-27oo, z<o. (1) 

TT(R 1 + R 2 ) 

Here, p is the bulk pressure at density p of the small spheres and 700 is the surface tension for small spheres of density 
p at a planar hard wall. For both quantities quasi-exact expressions are available 0,0! 
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In obtaining these equations we have set 

(5 = a = 1 (4) 

and we will do likewise in all following considerations. 

In recent years, progress has been made in calculating the depletion force in hard systems by other theoretical 
means, such as simulations 0,0, integral equations Q and density functional theory (DFT) 0. In a recent paper 
Henderson reviews this analysis of depletion forces in hard fluids and points to a serious discrepancy between 
Derjaguin's analysis on the one hand, and the various theoretical approaches/molecular dynamics simulations on the 
other hand. Also there are some features of the density functional result which do not fit with the simulations either. 
These discrepancies had been gone unnoticed partly due to the fact that comparisons were made between depletion 
potentials which add a fair amount of uncertainty to the simulation data since force curves with very few data points 
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FIG. 1: A comparison between existing results for the force between two colloidal particles in the depletion region for the size 
ratio a = 10. Shown are results from molecular dynamics 0, DFT [f|, the superposition approximation (using density profiles 
obtained as in and the Derjaguin limit for solvent densities: (a) p = 0.6 and (b) p — 0.7. 

had to be integrated Q . Also the limits inherent in the Derjaguin assumption have not been analyzed convincingly such 
that partial agreement respectively disagreement with Derjaguin's result has not been taken seriously. Henderson's 
analysis concentrated on the depletion force between hard walls and hard colloids but applies equally well to the force 
between two colloids. 

Before we analyze these discrepancies, we briefly present the strategies of the various approaches to obtain the 
depletion potential. Let us denote by p(r; xi,x 2 ) the density distribution of small spheres around two fixed hard 
spheres at positions xi and x 2 . Then the depletion force on one big sphere can be obtained by summing over all small 
spheres the force between a single small sphere and the big sphere. By symmetry, the force will be directed along the 
axis joining the centers of the two big spheres and due to the hard sphere interactions the volume integral reduces to 
an integral over the surface of one big sphere. Its magnitude (negative for attraction, positive for repulsion) is given 
by 

F a (z) = 2ir{R 1 +R 2 ) 2 J d(cos6>) cos6» p(r; 0,x 2 ) , (5) 
[|r| =R! + R 2 , x 2 = (R 2 + z, 0,0)] . 

In simulations, just this formula is used. The superposition approximation also uses this formula and additionally 
assumes that p(r; Xi,x 2 ) can be obtained by superimposing the two density distributions p(r — x^ around one fixed 
hard sphere centered at xi and x 2 , respectively: 

p(r; xi,x 2 ) = - p(r-xi)p(r-x 2 ) . (6) 
P 

The density distribution around one big sphere could be determined by e.g. integral equation methods or by min- 
imizing a density functional. When presenting superposition approximation results, we will use DFT results using 
the Rosenfeld functional as these are of superior quality. The DFT method of Ref. |5j also arrives at the depletion 
potential (whose derivative gives the depletion force) by just using the density distribution around one single big 
sphere but circumvents the crude approximation, eq. ©, by making use of the potential distribution theorem (also 
known as Widom's insertion trick). The method (which we call insertion route DFT) is explained in Add. IbI On 
the other hand, the distribution p(r; xi,x 2 ) could be obtained directly using DFT (in line with Ref. @ we call this 
brute force DFT). It is numerically involved and only two studies exist in the literature, both for size ratios smaller 
or equal a = 5 jg,0]- Error bars on the results of Ref. are much too large to arrive at a sensible conclusion. The 
much improved results of Ref. [7j indicate no significant deviation between the depletion potentials calculated using 
the insertion route and the brute force method, respectively. 

For densities p > 0.5 discrepancies between the above mentioned treatments and the simple Derjaguin formula 
become apparent as is illustrated in Fig.^ For a — 10 and solvent sphere densities p = 0.6 and 0.7 we show MD data 
Q, Derjaguin's result, insertion route DFT data calculated as in Ref. p| and data obtained from the superposition 
approximation. The deviation from Derjaguin's straight line is most obvious near 2 = 1, i.e. near where just one small 



3 



sphere fits between the two large spheres. The MD results seem to follow a straight line with a slope smaller than the 
one in Derjaguin's expression, p, but with a characteristic rounding off near z = 1 which always overshoots the DFT 
data (see also Fig. 7 in Ref . |j| for another simulation) . The DFT results show a flattening off which is characteristic 
for p > 0.5 and a > 10. The same behavior is seen in results using bridge diagram corrected HNC integral equations 
0. Finally, the superposition approximation produces a straight line with Derjaguin's slope but with a big offset. 
Using density distributions p(r — R) from less precise methods (integral equations with Percus-Yevick and Rogers- 
Young closure) offset and slope of the straight line are changed considerably such that these results fitted the MD 
data quite well. This led the authors of [j| to the erroneous conclusion that the superposition approximation is quite 
succesful in predicting p(r; Xi,X2). From the present results, it is clear that the superposition approximation does 
not constitute a good model of the force in the depletion zone. 

At first glance one is inclined to blame the discrepancies on the finiteness of a. After all, Derjaguin's result is 
supposed to be valid for a — * oo. Here the first problem arises: regarding this limit, Henderson gives an argument 
(which will be critically examined below) that deviations to Derjaguin should only occur for z > 1 — l/(4i?i + AR2) 
(z > 0.955 for a — 10). This is clearly not the case for all results as can be seen in Fig. ^ Moreover the insertion 
route DFT results do not converge to the Derjaguin limit for higher a 0. The second problem lies in the fact that 
Rosenfeld's (or related) DFT usually gives density distributions around fixed objects (wall |9j, big spheres [10(, wedge 
[Tlj ') of such a high quality that they seem to parametrize MC/MD data also for higher densities. However, in the 
present case systematic discrepancies between the MD and the DFT results occur. A tentative first explanation why 
this happens lies in the possibility that insertion route and brute force DFT give substantially different results for 
a > 10 (remember, there is no apparent difference for a = 5 |7j). A second possibility is that the higher-order 
correlations which are captured only approximately by any DFT model become more and more important. In fact, 
we will present below a picture for the depletion force which reveals quite subtle packing effects between the colloids 
which emerge for larger values of a. 

There is an interesting consequence from all of this. Defining the depletion potential by 



we note that in the Derjaguin approximation W a (0) < W a (l) only for p < 0.68. Above that density, contact between 
the two big spheres is only a metastable minimum separated by a rather high potential barrier from the overall 
minimum which will be close to z = 1. So, for higher densities the colloidal particles would not stick to each other. 
Although according to insertion route DFT Wq, (0) — Wa(l) also increases with increasing p > 0.7, this quantity never 
changes its negative sign for physical densities 26] . 

Therefore we can formulate our questions: Does the Derjaguin limit already set in for a « 10? If not, why? What 
is the source of discrepancy between DFT and MD/MC? As Rosenfeld's DFT is now being used in studies of solvation 
forces for colloidal particles in liquids with interactions other than hard sphere |l2j, the understanding of its limits 
for hard spheres is crucial. 

The remainder of the paper is organized as follows. In order to have a self-contained presentation, the Derjaguin 
limit for the depletion force is derived via (i) a force and (ii) energy analysis and (in) exact relations from statistical 
mechanics. This section contains nothing new and leans heavily on the presentation in Ref. To shed light on the 
onset of the Derjaguin limit, we will rederive it in a slightly different way in the following section and thus show that 
it is not valid when the colloids are separated by z ~ 1. This will define the new annular slit approximation. The 
depletion regime z < 1 is then analyzed in terms of an effective two-dimensional system of small disks which builts 
up in the annular wedge between the colloids. Using scaled particle theory in two dimensions, we derive an expression 
for F a (z) which for a — > 00 recovers the Derjaguin expression, although at a smaller rate as Henderson anticipated. 

For a = 10, the results of this analysis point to a flaw in the insertion route DFT treatment and show better 
agreement with the MD data. Therefore we will examine the insertion route DFT results closer and improve upon 
them by imposing test particle consistency (see App. lBjl . The equations obtained can also be viewed as RHNC integral 
equations with the bridge diagrams calculated from the hard-sphere density functional. Therefore, results from any 
integral equation closure can be viewed as being akin to insertion route DFT calculations. The quality of the density 
functional is then closely related to the quality of the bridge function approximation. 

Using test particle consistent DFT, we find no convergence to Derjaguin's result for a up to 100 and an increasing 
difference to the annular slit approximation (which becomes more reliable for increasing p and a). With the premise 
that the Dejaguin limit is reached in a non-singular way, we arrive therefore at the conclusion that insertion route 
DFT (and likewise integral equations) are unreliable for a > 10 since they miss some of the packing effects of the 
small spheres between the large colloids. 

In the last section we briefly comment on the possibility that non-analytic contributions might prevent a smooth 
transition to the Derjaguin limit. 




(7) 
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II. DERJAGUIN APPROXIMATION 



A. Force analysis 

The geometrical arrangement of the two colloids is shown in Fig. The boundary of the exclusion zone for the 
centers of the small particles is indicated by the dashed lines, thus the exclusion zone are two (possibly overlapping) 
spheres of radius R' — Ri + R^- The depletion force between the two large spheres is obtained by summing local 
pressures over the area of one (exclusion) sphere, 



F a (z) = 2nR' 2 J 



d{cos9)p loc (6) . (8) 



The Derjaguin approximation consists in replacing the local pressure by the solvation force per unit area of a planar 
slit with width I where the width refers to the minimal distance between the excluded volumes of the walls: 

Ploc(0) « /oo(0 • (9) 

Here I is the horizontal distance between the two (exclusion) spheres corresponding to the angle 8, see Fig. For 
very large a, this approximation is certainly justified, since locally the geometry resembles the planar slit. Using 

I = 2R' - 2x - (1 - z) , (10) 
cos 9 = x/R' , (11) 

we transform eq. iJSJ into 

/•OO 

F a (z)=nR' d//oo(0 (12) 

Jz-l 

The upper limit in this integral over the slit width has been replaced by infinity since it can be assumed that the 
solvation force approaches its limiting value foc(l — * oo) = when I is just a few a, which should be considerably 
smaller than R' . Now the solvation force per unit area is defined by 

M0 = -^, d3) 

where 7(0 is the excess grand potential (i.e. bulk grand potential subtracted) of the system of the two parallel walls 
which define the slit. Using the fact that 7(00) = 2700, where 700 is the surface tension of a single hard wall in a sea 
of small spheres we find 

F a {z)=nR'( 1 (z-l)-2 loo ) . (14) 

If z < 1 (i.e. no single small sphere fits into the slit) the surface tension arises from the release of free volume to the 
small spheres, 7 = p(z — 1). Thus 

F a {z) = 7r#(p(z - 1) - 2 7oo ) , z<l. (15) 

The depletion force is seen to depend only on the (hard sphere) pressure p and surface tension 700 for which we possess 
accurate approximations, see eqs. (J2GJ). 

B. Energy analysis 

Following Henderson |8| we can arrive at the Derjaguin approximation also by an analysis of the grand potential 
which can be decomposed into a 'volume', a 'surface area' and a 'line' term according to 

/•1 

?'2 



Cl(z) = -pV s (z) + 2ttR' 2 J d(cos 0)7(0 (16) 

P OO 

= - P V s (z)+ loc A s (z)+nR' <fl( 7 (i) - 2 7oo ) . (17) 








FIG. 2: (a) Geometrical definitions for two colloids in the depletion regime. The dashed lines indicate the surfaces of the 
exclusion spheres. Their overlap defines a volume Vg and an overlap surface area A' a . (b) Modelling the annular wedge part 
with widths I < lo by an annular slit of width e in which the solvent gas is effectively two-dimensional. 



Here, V s is the volume available to the small spheres (i.e. outside the two (possibly overlapping) exclusion spheres) 
and A s is the corresponding surface area of the two (possibly overlapping) exclusion spheres, 



v.(z) =v - v: = V - ~7ri*' 3 + \{z - 1) 2 (r> + L-Tj 



(18) 



A s (z) = A - A' s = 8ttR' 2 + 2irR'{z- 1) . (19) 

(Vq is the total system volume) and the last term is the 'line tension' contribution independent of the length of the 
overlap circle 27ryo and thus independent of z. Using 

dfl 

F a (z) = -— 20 
oz 

and retaining only the leading terms in 1/R', one arrives at the Derjaguin result, eq. Ijl5|l . 



C. Statistical mechanical analysis 

Summarizing Henderson's analysis, let us consider a fixed big sphere surrounded by a bath of small spheres. This 
sphere exerts an external, hard-body potential Vf Kt on small spheres [i — 1] and on any other large spheres [i = 2]). 

The work to insert another big sphere at distance z from the first sphere is given by —c^iz). Here, <4 (#) is the 
one-body correlation function given by 

4 1} (z) = log(p 2 (z)A 3 )-/i 2 (z), (21) 

A2O) = M-V^{z). (22) 

In these expressions, A 2 is the de-Broglie-wavelength of the big spheres and /x 2 is their chemical potential. For our 
configuration of interest, p 2 (z) is the one-body density profile of big spheres on another big sphere in the bath of 
small spheres in the dilute limit (/i 2 — * — 00). 

The first two equations in a hierarchy of functional derivatives of the grand potential are 

-ft(x) , (23) 



<5/2i( x ) 
J/I i (x)%(y) 



p i (x)p i (y)(^(x,y)-l)+ (0i (x)5( 3 )(x-y)^ . (24) 



6 



Let us assume that the center of the first, fixed sphere defines the origin of the coordinate system, and the coordinates 
of the center of the second big sphere are given by x = {2R 2 + z, 0, 0}. Since the external potential V2 vanishes for 
z > we find, using the above equations and the definition of the depletion force, eq. ifSHji. 

F a {z)=c^'{z) = EM (z>0 ), (25) 

P2(Z) 

In the dilute limit {p 2 — y 0) this expression simplifies to 

d 3 y Pi(y)(gi2(x,y)-1). (27) 

By virtue of the derivative of the external potential Vj 6 *' (exerted by the fixed sphere on the small spheres) only the 
surface of the exclusion sphere (i.e. one of the dashed lines in Fig.EJ) contributes to the integral. Since pi(y) is the 
density profile of small spheres around the fixed sphere, and pi (y) gi 2 (x, y) = p(y; 0, x) is the density profile of small 
spheres around the fixed sphere but with the second sphere fixed at position x, we recover eq. JSJ: 



F a (z) = 2ttR 12 J d(cos9)cos9 (p({R', 9}; 0, x) - pi(R')) , 



(28) 



/COS # m i n 
d(cos 9) cos 9 p({R', 9}; 0, x) , (29) 

cos^l 1 -^ 2 ^ [ Z z <\] (30) 

Although the integral in eq. (JSJ extends over the whole surface of the big sphere, the contribution of surface elements 
with azimuthal angle 9 < # m i n is zero since the density vanishes there. 

To identify the Derjaguin limit it is useful to keep the second term in the brackets on the rhs of eq. (|28|l (writing 
for the contact density of small spheres at a single large sphere p s = pi(R'j): 

/COS # m i n /-COS # m i n 

a da {p({R', a}; 0,x) - p s ) + 2ttR' 2 p s J a da (31) 

fCOS 6> min / ?r _ 1 \ 

= 2irR /2 J ada {p({R' , a};0 lX ) - p s ) + nR' p s {z - 1) N + J . (32) 

Changing integration variables from a = cos 9 to / according to eqs. (|10ll 1|) and identifying the upper limit with 
infinity, we find 

F a (z<l) = irj^ dl(R'-(l + l-z))(p({R\l};0 7 x)-p s )+nR'p s (z-l)(^l + ^^j . (33) 
The "microscopic" Derjaguin approximation consists in setting 

cos 9 (p({R', I}; 0, x) - p s ) = (l - ' + 2 ^7* ) (p(i R '' l }' °' x ) ~ P') ~ ^»(0 ~ P» ' ( 34 ) 

where p w (l) is the contact density at one wall in a planar slit of width I and p w is the contact density at a single 
planar wall. From statistical mechanics we furthermore know [l3l Il4| 

p w (l)-p w - /°°(0 = -^ ) (35) 

Ps = p + ^r + ltR- (36) 

Putting the last two equations into eq. (|33|l we find the Derjaguin result as the leading order in R' ', and we can 
identify the finite-size correction of first order to it (the surface tension on the (exclusion) sphere with radius R' can 
be approximated in first order by the surface tension on a hard wall, 7^ rj 7oo): 

F a ( Z < 1) W ^crjag + nRl 2 7oc + V ^ _ y (3?) 
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Interestingly, the finite size corrections predict a smaller slope for the force curves (~ 10% for a — 10) and a slight 
deviation from linearity which affects the curve only for z — > 0. We note that the considerations of Ref. 15] (their 
wedge approximation) would modify our finite-size corrections by F a /(wR') — > F a /(irR') — "foo/R 1 , i.e. the slope 
corrections would be mitigated. In any case, the qualitative behaviour for z — ► 1 remains unchanged. 



III. ANNULAR SLIT APPROXIMATION AND DERJAGUIN LIMIT 



The microscopic Derjaguin approximation of eq. I|34|) asserts that - apart from the geometrical factor cos0 - all 
annular wedges that are formed between the two large spheres for z < 1 are equivalent, i.e. the contact value of 
the density on the spheres can be described by a single function, namely p w (l)- At first glance, there is a physical 
difference between these wedges: At z = 1 the spheres on one "side" of the annular wedge can scatter with the spheres 
on the other side , as opposed to smaller values of z. Henderson argues that for small values of I, an effectively 
two-dimensional ideal gas of small spheres forms between the two colloids since the limiting three-dimensional density 
Pw{l — > 0) stays finite and therefore an effective 2D density p 2 d ~ I Pw vanishes. Therefore, scattering from one side 
of the wedge to the other should be negligible unless zero separation between the colloids occurs for radial distances 
yo < 1/2. Since ^Ri(l- z)R', it follows that the Derjaguin approximation should be valid for z < 1 — 1/(4R'). 

Is the concept of a nearly ideal 2D gas really valid in the annular wedge? For narrow slits with finite I we consider 
the small I <C a expansion of p w (l) [ljj: 

fhu{l) = — 7 ■ (38) 

P 1 exp( — p cx ) — 7T I 

Here, /i cx is the excess chemical potential of the small spheres and p the corresponding density. However, this limiting 
behavior is valid only for very small I. Rather one should study the effective 2D density in a slit defined by 

P2d(0= [ dl'p(l')=p l + T(l) ) (39) 
Jo 

with r(Z) defining the coverage in the slit. A Derjaguin-like estimate for the average 2D density in the annular wedge 
up to a maximal parallel distance Iq of the exclusion spheres follows: 

Pf A (h) = y f° P?J(l) = y f° T(l)dl + \p*h • (40) 

'0 Jo '0 Jo 1 

We can gain access to this quantity by using DFT again. As explained earlier, minimizing the Rosenfeld functional 
in the presence of an external field gives rather accurate density distributions. Carrying out the minimization in the 
presence of the two hard walls which define the slit gives us the explicit density distribution in the slit from which 
the surface tension 7, the coverage and the average 2D density as functions of I can be calculated. The results for two 
medium densities are shown in Fig. [3J It is seen that the coverage and therefore also the average 2D density quickly 
reaches the level of 2 T ^ (twice the coverage on a single planar wall) and therefore the 2D gas between the spheres is 
far from ideal. We also notice that the surface tension falls quickly to 2700 and then shows moderate oscillations 
around that value. 

Now, in order to formuate an alternative derivation of the Derjaguin limit we replace the last part of the annular 
wedge with I < Iq by an annular slit of width 1+e where the spheres can only move perpendicular to the z-axis, see 
the right panel of Fig. [21 The spheres in the slit can then be viewed as a system of hard disks. The surface tension 
in this fictitious slit is 7(^0) and its surface grand potential is written as 

fijrai = 7('o) Avcd + a(yo) 27n/o + ■ • ■ ■ (41) 

where A wc d is the one-sided area of the wedge and we have introduced a line tension term o~(yo) which describes the 
interaction at the inner boundary of the wedge. Now for large R' we have 



yo m y/R'(l - z) , A wed « A Q { yi )+nR'{z-l) , (42) 
where the area of a spherical cap is A Tty\. Now the depletion force has three contributions: 



F a (z<l) = (P+^) ^-^L+ I foo{ i )d i : N:!) 




FIG. 3: Surface tension y(l) and the negative of the coverage T(l) for planar slits of width I and two solvent densities: (a) 
p — 0.5 and (b) p — 0.7, as obtained from minimizing the Rosenfeld functional. The dot-dashed line shows the surface tension 
7(7 — > oo) = 2700. The (negative of the) average density refers to the Derjaguin like approximation of the 2D density in the 
annular wedge, eq. 14011 . 



which arise since we have split the original Derjaguin integral, eq. (|12|) . according to J°^_ 1 ■ ■ ■ — J^_ 1 ■ ■ ■ + L° ■ ■ ■+ ■ ■ ■ 
and have incoroprated the finite-size correction of eq. (|37|l . Simplifying eq. 143(1 using the geometrical relations in 
eq. we find, 



F a (z < 1) = ttR' 



27o 



R' 



4R' 



•7(io) + ^+</(!/o)-(27oo-7(io)) 



m 



(44) 



We have recovered the finite-size corrected Derjaguin result, eq. I|12fl . plus some line tension contribution which allows 
us to formulate Henderson's hypothesis of the equivalence of the annular wedges as 



lim = l im a '( yo ) = 

3/o-»0 y yo->0 



(45) 



This is not trivial at all. Rather, if cr(yo) or a ' iuo) g° to a finite value as yo — > we would expect the Derjaguin limit 
to fail for z — ► 1. 



A. The effective 2D gas in the annular slit 

We can gain access to the line tension function a(yo) by exploiting the nature of the quasi-2D gas in the wedge. 
Since the small spheres can only move in the plane perpendicular to the axis joining the centers of the two colloids, 
we rewrite the surface grand potential of our fictitious slit as 

O sur -> fi 2 d = -J32d(P2d)^wcd + 72d(P2d; 2/0 ) 27TJ/0 , (46) 

appropriate for the 2D "volume" and "area" grand potential contributions of a system of hard disks. Now let us 
think about what is physically happening when the second colloid approaches the first one at distances z < 1: A 
circular cavity forms in the center of the quasi-2D gas which cannot be reached by the centers of the solvent spheres. 
Therefore we can express the last equation as 

^2d = -P2d A {yi) + Mcav(yo) ( 47 ) 

Here, the first term is the ^-independent 2D "volume" term (Ao(yi) « ^y\, see the right panel of Fig. |2J), and we 
have introduced fJ- C nv(yb), the work needed to create a cavity of radius y' . From the reasoning above we would expect 
that y'o — yo- However, in our calculations of surface tensions and coverages in slits (see Fig.^J) we have seen that for 
very small slit widths (I < 5) p^ A — > 0. The limiting distance 5 can be estimated from the small I expansion of the 
contact density in slits, eq. (|38|) : 

5( P )«p- 1 exp(- / / x ) . (48) 
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This is indeed a small length compared to a: (5(0.5) « 4 ■ 10 2 , (5(0.7) 
wedge up to distances S must be added to the cavity, and therefore 



Vo 



y/R'((l-z) + 6) 



1 • 10~ 3 . But the depleted area in the annular 



(49) 



The problem of the insertion energy of an additional cavity was the starting point of scaled particle theory; here 
we can use the two-dimensional version Il7l to obtain 



Mcav(?/o) 
P2d 

72d 
£2d 



P2d T^v'i + 72d 2nyb + £ 2d 
-log(l -np 2d y' a 2 ) 



(y' > 1/2) 
(2/o < 1/2) 



V2d 



7T (1 - Tl2d) 2 

2 nld 



71" (1 - »?2d) 2 
1 - 277 2 d 



"?72d 



(1 - %d) 



l0g(l - ?72d) 



(50) 
(51) 
(52) 
(53) 
(54) 



The contact to the original surface energy of the slit, eq. I|41|) . is made by setting 7(^0 ) == — P2d, thus it follows that 
ciya) ~> 72d- To obtain numbers, we simply choose Iq such that 7(^0) = 2700. Using eq. for the 3D surface tension, 
we can determine p 2 d as a function of p, the 3D density of small spheres. Remember that physically it would be also 
quite sensible to identify the 2D density via eq. H39fl . P2d = P2d- A quick glance at Fig. [21 assures us that the two 
definitions of p 2 d are quite consistent with each other |27|. 

One might be concerned that the validity of eq. I|50|l is limited for intermediate disk sizes y' Q = 1 ... 2 (where we 
would need it if calculating numbers for a — 10, say) as scaled particle theory is by construction only an interpolation 
between the known analytical behavior of A t cav(j/o) l° r Vo < 1/2 011 the one side and y' — > 00 on the other side. In 
view of lack of appropriate data in the literature we have performed a quick MC check for two densities p 2 d = 0.4 and 
0.6, the results are shown in appendix ^] From these results it follows that scaled particle theory is precise enough 
for our purposes. 

The final result for the depletion force, following from eqs. (|43I44|) and the considerations of the previous paragraphs, 
takes the form: 



F a (z < 1) 1 



ttR' 



27oo \ dV s 1 dp 

C 



ttR' \ P+ R' ) dz irR' dz ' 



270 



R 1 



1 



P2d 



l-z + 6 

2R' 

1 - np2dy'o 



1-z + S 
2R 7 



(Vo < 1/2) 



(Vo > I/ 2 ) 



{1-z + S) (R' 



1-z + S 



(55) 
(56) 

(57) 
(58) 



Here, for the sake of completeness, the exact geometrical expression for the cavity radius y' is given. The expression 
in eq. I|49|l is the leading term in an expansion of y' with respect to R' . 
Most remarkably, it follows for the force at z = 1 (to first order in 1/R'), 



F a (z = l) 
irR' 




(SR 1 < 1/4) 
(8R' > 1/4) 



(59) 



i n 





FIG. 4: The scaled force between two colloidal particles in the depletion region z < 1 for a size ratio a = 10 and for three 
solvent densities: (a) p = 0.5, (b) p = 0.6 and (c) p — 0.7. Comparison between the annular slit approximation (eq. I56H ). MD 
data and test particle consistent DFT based on the White Bear (WB) functional. 
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0.4 


0.22 


0.18 




0.28 


0.22 


0.22 


0.31 


0.5 


0.34 


0.04 




0.44 


0.41 


0.36 


0.62 


0.6 


0.46 


8 • 10" 


3 


0.49 


0.57 


0.51 


1.15 


0.7 


0.59 


1 ■ 10" 


3 


0.60 


0.74 


0.62 


2.05 


0.8 


0.71 


5 ■ 10" 


5 


0.71 




0.64 


3.55 


0.9 


0.81 


8 • 10" 


7 


0.81 




0.62 


6.01 



TABLE I: Results for the scaled depletion force at z = 1, F a /ljiR'), for the annular slit approximation, from MD simulations 
0, test particle consistent DFT and the Derjaguin approximation (a = 10). Note that the annular slit approximation predicts 
that the value of the scaled force at this point is essentially given by p2d, the effective 2D density. 
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IV. RESULTS FROM THE ANNULAR SLIT APPROXIMATION AND COMPARISON WITH DFT 

AND MD 

Having obtained a closed expression for the force in the depletion region, eq. Q56I1. we can compare results to the 
available MD data and to the DFT results. Instead of using the results from Ref. J5| we apply the bridge functional 
formalism of Ref. 18] to obtain results which are "one test particle consistent". In effect, the equations for the 
depletion potential are transformed into RHNC type equations where Rosenfeld's density functional (or extensions 
thereof) is the generating source for the bridge diagrams. For more details we refer to App.[B] where we have outlined 
the procedure and compared to the previous DFT results. Summarizing the results from the appendix, the self- 
consistency for one test particle gives a relatively small shift of the depletion force for z < 0.6 . . . 0.7 which is always 
upwards. This adds up to a 10 ... 20 % correction upwards for the depletion potential at contact, W a (z = 0). For 
z > 0.7, the results are quantitatively almost unchanged. Especially the failure of the previous results to converge to 
the Derjaguin limit remains unaltered. 

For the largest ratio a = 10 where simulation data are available, we show results for the depletion force in Fig. 0] 
Values for the force at z = 1 are compared in Tab. [i] In general, the agreement between the simulations and our 
annular slit approximation is surprisingly good. The approximation follows the trend of the simulation data to produce 
a maximum in the depletion force for z < 1 and p > 0.5. A pronounced maximum is absent in the DFT results for 
the densities shown. Despite the better agreement of the annular slit approximation with the simulation data, it 
is hard to tell whether for this size ratio a there is already a serious problem with DFT. First, there are no error 
bar estimates for the MD data available, and secondly, the approximation suffers from possible errors due to a finite 
number of particles in our idealized annular slit. This number can be estimated by 

N s = p 2 dAvcd ~ P2d^R'k (60) 

Indeed, since Iq rs 0.2 . . . 0.4 (see Fig. |2J) N s < 5 for all densities! According to this estimate, our considerations 
should become increasingly reliable for larger /?2d and larger a. Better correspondence with the simulation results 
with increasing p2d is indeed observed, see Fig. 

For larger a, the discrepancy between the annular slit approximation and the DFT results becomes striking. We 
show this in Fig. [5] for two size ratios, a — 10 and a = 100. As we have explained, the annular slit approximation 
can be expected to become more accurate for larger a and it has the correct limiting behavior, so the conclusion 
would be that DFT becomes increasingly unreliable for a > 10. Although not shown in the figure, there is already 
a substantial difference for a ~ 20, say. Thus one should regard with extreme caution the claim of |f| that insertion 
route DFT can be expected to be rather accurate also for size ratios larger than 10. A similar claim made about a 
bridge diagram improved HNC treatment of the depletion potential (see Ref. Q) should also be treated with caution 
as the HNC results show similar defects as the DFT results. Recall that the improved, test particle cosistent DFT 
results shown here can be viewed as HNC results with bridge diagram corrections supplied by the density functional 
and both methods can be formulated in the language of insertion route DFT. 

Finally we calculate the quantity W a (Q) — W a (l) (for a — 10) which is roughly the depletion potential difference 
between colloid contact and the first minimum for medium to high densities. The results are collected in Tab. ITU We 
see that the previous DFT results predict that the potential at colloid contact is minimal for all values of p. This 
finding is not changed by imposing test particle consistency; only the absolute value of the potential difference is 
reduced somewhat. The annular slit approximation predicts that the absolute minimum jumps to z ~ 1 for p ss 0.83, 
still far away from the Derjaguin value 0.68. For a = 100, however, the jump of the absolute minimum occurs at a 
density of 0.73, according to the annular slit approximation. 

V. SUMMARY AND CONCLUSIONS 

In this paper, we have analyzed the depletion force between hard colloids in a solvent consisting of small hard spheres. 
We started from an already previously observed disagreement between simulation data/results from Rosenfeld's DFT 
and the Derjaguin limit. Albeit first derived on purely phenomenological grounds, equilibrium statistical mechanics 
strongly supports the validity of the Derjaguin limit for large size ratios a between colloids and solvent spheres. 
The disagreement between simulation data/DFT and the Derjaguin limit in the depletion zone for a = 10 could be 
explained by an effective approximation which analyzes the structure of the solvent between the two colloids in terms 
of a (fairly) dense 2D gas of hard disks. The depletion force near the onset of the depletion zone (i.e. where exactly 
one solvent sphere fits between the colloids) is mainly determined by the force to create a disk cavity in the effective 
2D gas. Requiring the Derjaguin limit for a — > oo, there are no free parameters for the 2D gas. The agreement with 
simulation data is very good even at the relatively small size ratio a — 10. 
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FIG. 5: Scaled depletion force for size ratios a = 10 and a — 100: Comparison between the annular slit approximation (gray 
curves) and test particle improved DFT based on the White Bear (WB) functional. The solvent density is p — 0.6. Note that 
the annular slit approximation approaches the Derjaguin limit quite slowly, nevertheless the Derjaguin limit is not reached at 
all by the DFT results. 



p W a (0)-W a (l) 

test particle annular slit 
DFT \5] improved DFT approximation Derjaguin 



0.4 -3.65 -3.50 -3.15 -3.30 

0.5 -4.32 -4.01 -3.46 -3.35 

0.6 -4.69 -4.16 -3.50 -2.27 

0.7 -4.70 -3.84 -2.89 0.87 

0.8 -4.39 -3.19 -1.02 7.79 

0.9 -3.77 -2.38 3.12 21.64 



TABLE II: Depletion potential difference W a (0) — W a (l) = F a (z)dz (a = 10) calculated using insertion route DFT Q, test 
particle improved DFT based on the White Bear functional, the annular slit approximation and the Derjaguin approximation. 



For higher size ratios a > 10 no simulation data are available, and existing DFT results showed no convergence 
towards the Derjaguin limit. Imposing test particle consistency we calculated improved DFT results which however 
did not alter their large a behavior. Already for a — 20 the disagreement of the DFT depletion force with the results 
of the effective annular slit approximation becomes pronounced. We conclude that the limit of reliability of insertion 
route DFT is reached for a = 10. Through the test particle consistent calculations we have shown that insertion 
route DFT and integral equation approaches are methodologically equivalent to each other; the single variants differ 
in their choice for the bridge functional. Similar limits for the reliablility can therefore be also expected for integral 
equations. Likewise a similar limit will apply if one treats solvophobic colloids in attractive fluids with repulsive cores 
using Rosenfeld's DFT for hard sphere reference systems. To reach larger size ratios, the analysis of the annular slit 
approximation could also be extended to this case. 

We have shown that quite subtle packing effects between the colloids play a role in determining the depletion force 
for medium to large size ratios. In light of this it is actually amazing that insertion route DFT (which explicitly needs 
only the density distribution around one colloid) captures most of the effects and only misses the intricate effect of the 
quasi-2D gas. We emphasize again that the insertion procedure is formally exact but we possess only approximate 
expressions for the hard sphere density functional whose functional derivative is needed for the insertion procedure 
to work. The two variants investigated herein, the Rosenfeld and the White Bear functional, are - despite being very 
precise - not exact. One deficiency, if not the main, lies in the bulk direct correlation functions of order n > 2: 
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they are zero outside the hard core according to the functionals, but from simulations and integral equations we know 
otherwise. Requiring test particle consistency has fixed this shortcoming for c^ 2 \ but for all higher order (n > 2) 
correlation functions one should require consistency for n — 1 fixed test particles. Part of the problem for n = 3 is 
thus the determination of the density profile with the two colloids fixed. Therefore, it would be interesting to see how 
brute force DFT fares for size ratios a > 10, i.e. whether one could observe the oscillatory packing in the annular 
wedge directly and how it is related to . 

Throughout the paper we have argued that the Derjaguin limit for the depletion force is meaningful. The annular 
slit approximation, eq. I|56|l . predicts the leading correction to the force oc a -1 / 2 . Interesting enough, this is a non- 
analytic term but the Derjaguin limit is still reached continuously in the variable 1/a. Since the depletion force is also 
connected to an integral over the surface densities, eq. (J5J, this constitutes a hint that the density profile also contains 
non-analytic contributions in 1/a. In fact, whereas there are good arguments that the density profile around hard 
convex objects should have an analytic expansion in terms of the curvatures , such an analyticity requirement 
does not hold for profiles around non-convex objects (such as the two-colloid configuration in the depletion region). 
Consequently, there is the possibility that for 1/a — ► the surface densities do not reach the Derjaguin limit of the 
surface densities in a planar slit configuration, eq. Ij38(l . Possible singular contributions can not be understood with 
the current theories due to entropic arguments. A finite difference between the surface density and its Derjaguin limit 
would also constitute a surprise since it would point to small sphere correlations which are much larger than the bulk 
correlation length and allude to a mysterious phase transition. Nevertheless the non-analyticity of density profiles 
around curved objects is an extremely interestig subject in itself which is currently being explored 20]. 
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APPENDIX A: QUALITY OF SCALED PARTICLE THEORY FOR HARD DISKS IN 2D 

To check the reliability of the 2D scaled particle theory predictions for the insertion energy of a cavity with radius 
2/o, Wcav(yo)j we performed Monte Carlo tests for two medium densities, p2d = 0.4 and p2d — 0.6. The insertion energy 
is given by 

Wcav(yo) = logpo(2/o) , (Al) 

where po(j/o) is te probability to find no center of a disk within a circle of radius j/o around an arbitrary point. 
Therefore we chose for each Monte Carlo move a new, random point, around which we checked the latter condition 
[2l| . The results, obtained for 4000 disks and roughly 10 s moves, are depicted in Fig. EI The MC results are compared 
to scaled particle theory, and it is found that its simple prediction 

Wcav(2/o) = 7T2/0 P2d + 27TJ/ 72d (A2) 

with 72d independent on the cavity radius is extremely good except in the vicinity of yo — 1/2. At this point, exact 
analysis demands that the third derivative has a singularity pq. 

^( Sb _l/2 + )«( tt) -l/2)- 1 /a. (A3) 

This leads to a square root like cusp in the second derivative as can be seen in the MC results. Scaled particle theory 
ignores this cusp which is not too bad an approximation since the cusp quickly relaxes to a constant which is the 
pressure of the disk system. Apart from this effect of the non-analyticity of w cav , scaled particle theory is sufficiently 
precise for our purposes. 



APPENDIX B: TEST-PARTICLE CONSISTENT DFT 



In this section we shortly explain the route taken by ; 5| for obtaining the depletion potential, and introduce correc- 
tions thereto by requiring test particle consistency. Numerical results for both methods are presented. Throughout 
this section (3 = a = 1. 




FIG. 6: Monte Carlo results for the insertion energy (and its first two derivatives) of an (exclusion) disk of radius yo into a 
system of hard disks with diameter 1. The thick black lines are the scaled particle predictions, the rugged gray lines are raw 
data. The disk densities are: (a) p2d = 0.4 and (b) p2d = 0.6. 



Suppose we have a mixture of two hard species with particle radii R\ and i?2 • Species 2 refers to the colloids and 
the mutual interaction potentials are given by Uij (r) . We split the density functional describing the mixture into an 
ideal gas and an excess part, 

r[{(H(r)}] = ^d[{^(r)}] + ^ ex [{pi(r)}] . (Bl) 
The hierarchy of direct correlation functions is defined as 

c ? t ( Xl ...x n ) = - -= , \ , j ; ■ (B2) 

Specifically for n = 1 one can introduce an excess chemical potential functional, 

/ iHx;{ /0i (r)}] = -cf ) (x), (B3) 

which reduces to the usual excess chemical potential if the densities are constant, /i° x [x; {pi(r) = const}] = Pi X {{pi,o})- 
Suppose we have an inhomogeneous situation where one colloid is fixed. The depletion potential is then defined as 
the difference between the work needed to put another colloid particle into the system at position x on the one hand 
and at infinity on the other hand. Using the potential distribution theorem 0] we find 

W a (x) = lim (mHx; {Mr)}} - ft"({fto})) , (B4) 

P2— »0 

where we have assumed that limix!—^ pi(x) = pi,o- Note that p^l- ■ ■] depends in the required limit only on the 
density distribution pi(r) of species 1 around the colloid, i.e. before the second colloid is inserted. With a given excess 
functional at hand, the depletion potential is found by obtaining p\ (r) through grand potential minimization, 

SF 

= __- /ii + y i (x - (B5) 
-log(pi(x)) = Mf[x;{pi(x'),0}]-/if({p 1 , o ,0})+u 12 (x) , (B6) 
and then inserting this solution into eq. (|B4J| . 

1. Test particle consistency 

The depletion potential is the negative potential of mean force between a pair of colloids at infinite dilution. The 
following relation is valid: 

W r a (|x|) = -Iog fl a2(|x|;{p 1)0j 0})-U22(|x|) , (B7) 
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where 322 is the colloid-colloid distribution function in the limit of vanishing colloid density. On the one hand, this 
distribution function can be determined via the depletion potential described in the manner above. On the other 
hand, the excess functional defines the second-order correlation function Cy through eq. I|B2(I (for n = 2) which in 
turn can be inverted to give gij using the Ornstein-Zernike relation 

M|r|)-i 2) (|r|) = 5>*.°'*** C 8 ) (M)' (B8) 
fc 

MM) = MM)-1, (B9) 

^*cg } (|r|) = Jd 3 v'h ik {\v'\)c^{\v-v'\). (BIO) 

In general, both routes will give different results for an approximated free-energy functional. 

To make both routes consistent with each other, we proceed as follows Consider the equation which relates 
the bridge function bij to the distribution and direct correlation function, 

MM) = -iogM|r|)-MM) + E<°M^*4?(l r l)- (bii) 

k 

The various closures of integral equations follow by specifying a model for the bridge function, e.g. bij{r) — for the 
HNC closure. The bridge function can be generated by a bridge functional which we define to be the functional which 
contains all contributions beyond second order in a density expansion of the exact free energy functional around fixed, 
constant bulk densities = pi(r) — Apj(r): 



F[{pi(M)}] = ^id[{p fc (|r|)}]+^ ex (W)+Mf(W) Jd 3 rA Pi (v)- (B12) 

~ Jd 3 v Jd 3 v' 4 2) (r, r'; bulk)A Pi (r)Ap,(r') + j£[{p*(r)}] . (B13) 

Doubly occuring indices are summed over. To verify that the such introduced F~ indeed generates the bridge 
functions, we minimize the grand potential according to eq. i|B5() in the presence of the interparticle potential, 
Vi = Uij . 

= log P M. _ Ui ( r ) + $ * Apk (r) . (B14) 

OPi( r ) Pi,0 

Since pi ; o <7ij(|r|) = Pi(|r|) and pi t o hij(\r\) — Api(|r|) we have recovered eq. IjBllll upon the identification 



M r ) 



8T] 



br 



8 P i(r) 



(B15) 

Vi=Uij 



Up to now everything was exact but in order to specify a closure explicitly we assert that the true bridge functional 
can be approximated by the bridge functional of a reference model for which we possess an explicit form of JF CX : 

^[{PfcW}] * -^ f [{Pfc(r)}] - AC rcf ({PM}) Jd 3 r A Pl (r) + \ Jd 3 r Jd 3 r> cf ' rcf (r, r'; bulk)A^(r) A Pj (r') .(B16) 

A remark is in order here. Note that eqs. (|B12IB16fl together define a new functional which is now test particle 
consistent, i.e. the inversion of the Ornstein-Zernike relation gives the same result as an explicit determination of the 
distribution functions through the density profiles around test particles. This consistency holds regardless of the form 
of the interparticle potential and of how good or bad the choice of the reference system is. Of course, the reference 
system of choice is again hard spheres described by Rosenfeld's functional or a recently improved version, the White 
Bear functional 

In the limit P2 — > 0, relevant for the determination of the depletion potential, the equations for the distribution 
functions decouple. For <?n we have to solve the following equations (r = |r|): 

hu(r) - 4? (r) = p h0 h n * c$ (r) , (B17) 
-logsu(r)-un(r) - p 1)0 fcu * - $)(r) + (^[r-, {pi, g u (r), 0}] - ^({Pi.o.O})) . (B18) 




FIG. 7: (a) Depletion potential and (b) depletion force between two colloids for a solvent density of p = 0.8 and a size ratio 
a — 10: Comparison between test particle inconsistent and consistent results, using Rosenfeld's (RF) and the White Bear 
(WB) functional respectively. 



As expected, the colloids decouple and we are left with the equations for the one-component system. Employing the 
White Bear functional, we obtain one-component distribution functions which fit the MC data even better than the 
standard Verlet parametrization [22l l23| . The input of /in, en is needed to solve the next two equations for gi 2 (the 
normalized density distribution around one colloid): 

hu(r) - cf 2 \r) = p lfi h n * cf 2 \r) , 
-log #12 (r) - u 12 (r) = pi, hi2 * (4i ),re£ - 4f)(r) + (/^ x,ref [r; {p lj0 gu(r),0}] - fj, e *' Tei ({p lj0 , 0})) 

Having obtained h±2, C12, the depletion potential is simply given by 

W a (r) = -\ogg 22 (r) - u 22 (r) = p 1>0 h 12 * (c^ - c$)(r) + [r; { Pl . <?i 2 (r),0}] - p?' Iei ({p lj0 , 0})) (B21) 

Comparing the expressions for the depletion potential for test particle inconsistent and consistent DFT, 
eqs. (IB4IB21|) . wc sec that the main difference is buried in the first term on the rhs of eq. i|B21(l since p\ (r)|i ncons i stG nt ~ 

Pl,0 .912 Ml consistent • 

Results for the size ratio a = 10 reveal no huge differences between the test particle consistent and inconsistent 
calculations. For distances between the colloids z < 0.6 . . . 0.7 the consistent results give a somewhat higher force than 
the inconsistent results which adds up to a noticeable upward shift in W a (z = 0). Apart from that the differences are 
minimal, even at higher densities. For p^o = 0.8, we show the depletion potential and force in Fig. calculated with 
the Rosenfeld and the White Bear functional. 

Due to the accuracy in the one-component case, this choice of reference system has been extended to binary soft 
systems 0. Again, the agreement with simulation data is extremely good but the size ratio in the binary systems 
was well below 10. Only recently the depletion potential between soft colloids in soft fluids has been calculated for 
size ratios of about 10 using this method l- n general, test particle consistent DFT fares much better than any 
other theoretical method compared to the simulation data. Nevertheless, in the case of hard colloids in a LJ fluid 
discrepancies to the simulation data occur which show the same footprints as the deviations we observe here in the 
case of hard colloids in hard fluids. 
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